For execution on a local, multicore CPU with excess RAM we recommend calling
options(mc.cores = parallel::detectCores()).
To avoid recompilation of unchanged Stan programs, we recommend calling
rstan_options(auto_write = TRUE)
For within-chain threading using `reduce_sum()` or `map_rect()` Stan functions,
change `threads_per_chain` option:
rstan_options(threads_per_chain = 1)
Loading required package: gdata
Attaching package: 'gdata'
The following object is masked from 'package:stats':
nobs
The following object is masked from 'package:utils':
object.size
The following object is masked from 'package:base':
startsWith
Loading required package: bayesplot
This is bayesplot version 1.11.1
- Online documentation and vignettes at mc-stan.org/bayesplot
- bayesplot theme set to bayesplot::theme_default()
* Does _not_ affect other ggplot2 plots
* See ?bayesplot_theme_set for details on theme setting
Loading required package: stringr
Loading required package: dplyr
Attaching package: 'dplyr'
The following objects are masked from 'package:gdata':
combine, first, last, starts_with
The following objects are masked from 'package:stats':
filter, lag
The following objects are masked from 'package:base':
intersect, setdiff, setequal, union
Loading required package: ggplot2
Loading required package: loo
This is loo version 2.7.0
- Online documentation and vignettes at mc-stan.org/loo
- As of v2.0.0 loo defaults to 1 core but we recommend using as many as possible. Use the 'cores' argument or set options(mc.cores = NUM_CORES) for an entire session.
Attaching package: 'loo'
The following object is masked from 'package:rstan':
loo
Loading required package: hBayesDM
Loading required package: Rcpp
This is hBayesDM version 1.2.1
Attaching package: 'hBayesDM'
The following object is masked from 'package:bayesplot':
rhat
Loading required package: tidyr
Attaching package: 'tidyr'
The following object is masked from 'package:gdata':
starts_with
The following object is masked from 'package:rstan':
extract
# extract log likelihood for each choicelog_likelihood_main <-extract_log_lik(fit_main, parameter_name ="log_lik", merge_chains = T)log_likelihood_main_DU <-extract_log_lik(fit_main_DU, parameter_name ="log_lik", merge_chains = T)log_likelihood_PH_withC <-extract_log_lik(fit_PH_withC, parameter_name ="log_lik", merge_chains = T)#log_likelihood_PH_withC_DU <- extract_log_lik(fit_PH_withC_DU, parameter_name = "log_lik", merge_chains = T)# create logical vector coding if a column should be included in log_likelihoodx <-logical(stan_data$T) # create vector containing total number of trials * FALSEa <-c(stan_data$N+1:(stan_data$N*49)) x[a] <-TRUE# set logical value to TRUE if index in a# exclude 1st trial per subject as likelihood is identical and hinders loo estimation (causes Inf pareto k values)log_likelihood_main <- log_likelihood_main[,x] log_likelihood_main_DU <- log_likelihood_main_DU[,x]log_likelihood_PH_withC <- log_likelihood_PH_withC[,x] # exclude missing trialslog_likelihood_main <- log_likelihood_main[,log_likelihood_main[1,]!=-999]log_likelihood_main_DU <- log_likelihood_main_DU[,log_likelihood_main_DU[1,]!=-999]log_likelihood_PH_withC <- log_likelihood_PH_withC[,log_likelihood_PH_withC[1,]!=-999]#log_likelihood_PH_withC_DU <- log_likelihood_PH_withC_DU[,log_likelihood_PH_withC_DU[1,]!=-999]
2.3 Fit per model
Code
loo_main <-loo(log_likelihood_main)
Warning: Can't fit generalized Pareto distribution because all tail values are
the same.
Warning: Can't fit generalized Pareto distribution because all tail values are
the same.
Warning: Can't fit generalized Pareto distribution because all tail values are
the same.
Warning: Can't fit generalized Pareto distribution because all tail values are
the same.
Warning: Can't fit generalized Pareto distribution because all tail values are
the same.
Warning: Can't fit generalized Pareto distribution because all tail values are
the same.
Warning: Can't fit generalized Pareto distribution because all tail values are
the same.
Warning: Can't fit generalized Pareto distribution because all tail values are
the same.
Warning: Can't fit generalized Pareto distribution because all tail values are
the same.
Warning: Can't fit generalized Pareto distribution because all tail values are
the same.
Warning: Can't fit generalized Pareto distribution because all tail values are
the same.
Warning: Can't fit generalized Pareto distribution because all tail values are
the same.
Warning: Can't fit generalized Pareto distribution because all tail values are
the same.
Warning: Can't fit generalized Pareto distribution because all tail values are
the same.
Warning: Can't fit generalized Pareto distribution because all tail values are
the same.
Warning: Can't fit generalized Pareto distribution because all tail values are
the same.
Warning: Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.
Code
print(loo_main)
Computed from 36000 by 5583 log-likelihood matrix.
Estimate SE
elpd_loo -1497.2 40.6
p_loo 134.2 12.7
looic 2994.5 81.3
------
MCSE of elpd_loo is NA.
MCSE and ESS estimates assume independent draws (r_eff=1).
Pareto k diagnostic values:
Count Pct. Min. ESS
(-Inf, 0.7] (good) 4935 88.4% 228
(0.7, 1] (bad) 425 7.6% <NA>
(1, Inf) (very bad) 223 4.0% <NA>
See help('pareto-k-diagnostic') for details.
Code
plot(loo_main)
Warning in plot.psis_loo(loo_main): 0.29% of Pareto k estimates are Inf/NA/NaN
and not plotted.
Code
if (any(pareto_k_values(loo_main) >0.7)) { loo_main_clean <-loo(log_likelihood_main, k_threshold =0.7)print(loo_main_clean)plot(loo_main_clean)}
Warning: Can't fit generalized Pareto distribution because all tail values are
the same.
Warning: Can't fit generalized Pareto distribution because all tail values are
the same.
Warning: Can't fit generalized Pareto distribution because all tail values are
the same.
Warning: Can't fit generalized Pareto distribution because all tail values are
the same.
Warning: Can't fit generalized Pareto distribution because all tail values are
the same.
Warning: Can't fit generalized Pareto distribution because all tail values are
the same.
Warning: Can't fit generalized Pareto distribution because all tail values are
the same.
Warning: Can't fit generalized Pareto distribution because all tail values are
the same.
Warning: Can't fit generalized Pareto distribution because all tail values are
the same.
Warning: Can't fit generalized Pareto distribution because all tail values are
the same.
Warning: Can't fit generalized Pareto distribution because all tail values are
the same.
Warning: Can't fit generalized Pareto distribution because all tail values are
the same.
Warning: Can't fit generalized Pareto distribution because all tail values are
the same.
Warning: Can't fit generalized Pareto distribution because all tail values are
the same.
Warning: Can't fit generalized Pareto distribution because all tail values are
the same.
Warning: Can't fit generalized Pareto distribution because all tail values are
the same.
Warning: Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.
Computed from 36000 by 5583 log-likelihood matrix.
Estimate SE
elpd_loo -1497.2 40.6
p_loo 134.2 12.7
looic 2994.5 81.3
------
MCSE of elpd_loo is NA.
MCSE and ESS estimates assume independent draws (r_eff=1).
Pareto k diagnostic values:
Count Pct. Min. ESS
(-Inf, 0.7] (good) 4935 88.4% 228
(0.7, 1] (bad) 425 7.6% <NA>
(1, Inf) (very bad) 223 4.0% <NA>
See help('pareto-k-diagnostic') for details.
Warning in plot.psis_loo(loo_main_clean): 0.29% of Pareto k estimates are
Inf/NA/NaN and not plotted.
Code
loo_main_DU <-loo(log_likelihood_main_DU)
Warning: Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.
Code
print(loo_main_DU)
Computed from 36000 by 5583 log-likelihood matrix.
Estimate SE
elpd_loo -1515.1 42.0
p_loo 122.9 13.2
looic 3030.2 83.9
------
MCSE of elpd_loo is NA.
MCSE and ESS estimates assume independent draws (r_eff=1).
Pareto k diagnostic values:
Count Pct. Min. ESS
(-Inf, 0.7] (good) 4847 86.8% 899
(0.7, 1] (bad) 340 6.1% <NA>
(1, Inf) (very bad) 396 7.1% <NA>
See help('pareto-k-diagnostic') for details.
# extract log likelihood for each choicelog_likelihood_PH_withC_hierarchical_group <-extract_log_lik(fit_PH_withC_hierarchical_group, parameter_name ="log_lik", merge_chains = T)# create logical vector coding if a column should be included in log_likelihoodx <-logical(stan_data_hierarchical$T) # create vector containing total number of trials * FALSEa <-c(stan_data_hierarchical$N+1:(stan_data_hierarchical$N*49)) b <-c(stan_data_hierarchical$N*50+stan_data_hierarchical$N+1:(stan_data_hierarchical$N*49))c <-c(a,b) # concatenate both number vectorsx[c] <-TRUE# set logical value to TRUE if index in c# exclude 1st trial per subject as likelihood is identical and hinders loo estimation (causes Inf pareto k values)log_likelihood_PH_withC_hierarchical_group <- log_likelihood_PH_withC_hierarchical_group[,x] # exclude missing trialslog_likelihood_PH_withC_hierarchical_group <- log_likelihood_PH_withC_hierarchical_group[,log_likelihood_PH_withC_hierarchical_group[1,]!=-999]
# Load log_pd for model with fixed effect for grouplog_pd_PH_withC_hierarchical_group <-readRDS(file.path(out_path, "log_pd_kfold_bandit2arm_delta_PH_withC_hierarchical_group_n58.rds"))# exclude 1st trial per subject as likelihood is identical and hinders loo estimation (causes Inf pareto k values)log_pd_PH_withC_hierarchical_group <- log_pd_PH_withC_hierarchical_group[,x] # exclude missing trialslog_pd_PH_withC_hierarchical_group <- log_pd_PH_withC_hierarchical_group[,log_pd_PH_withC_hierarchical_group[1,]!=-999]